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We discuss various properties of the ground 
state of a Bose-condensed dilute gas con- 
fined by an external potential. We devote 
particular attention to the role played by 
the interaction in determining the kinetic 
energy of the system and the aspect ratio 
of the velocity distribution. The structure of 
the wave function near the classical turn- 
ing point is discussed and the drawback of 
the Thomas-Fermi approximation is ex- 
plicitly pointed out. We consider also states 
with quantized vorticity and calculate the 



critical angular velocity for the production 
of vortices. The presence of vortex states 
is found to increases the stability of the 
condensate in the case of attractive inter- 
actions. 
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1. Introduction 



In this paper we discuss some relevant properties of 
the ground state of a Bose-condensed atomic gas con- 
fined by an external potential. Our starting point is the 
Gross-Pitaevskii equation which gives the proper 
Schrodinger equation for the order parameter of an in- 
homogeneous dilute Bose-condensed gas at zero tem- 
perature. Using this equation it is possible to discuss 
various ground state properties of the system: the form 



of the atomic cloud, the role of the interatomic potential, 
the velocity distribution, and so on. All these quantities 
are essential for the interpretation of the recent experi- 
ments on Bose-Einstein condensation in ultracold alkali 
atom gases [1-3]. An important feature is the profound 
difference between systems interacting with repulsive 
and attractive forces. In the latter case, in particular, the 
stationary solution given by the Gross-Pitaevskii equa- 



537 



Volume 101, Number 4, July-August 1996 

Journal of Research of the National Institute of Standards and Technology 



tion is of metastable type. If the number of atoms is too 
large such a solution becomes unstable and the system 
collapses. In this paper we will also discuss some rota- 
tional properties of the system, in particular the struc- 
ture of vortices and the critical angular frequency 
needed to generate a rotational instability. 

The Gross-Pitaevskii equation for the order parameter 
ijjir) = {ipQ')) has the well known form [4]: 



2m 



V^ + VUr) + 



4'nA^a 



\^(r) 



ijjir) = fjujjir), (1) 



where Kxt is the external confining potential, which is 
usually chosen in the form of an anisotropic harmonic 
well. The role of interactions is accounted for by the 
non-linear term and is parametrized by the s -wave scat- 
tering length a . The quantity /^ is the chemical potential 
and is fixed by imposing the proper normalization, 
N=J p6r, to the density of the system p= |i/^|^. The 
Gross-Pitaevskii equation ignores interaction effects due 
to the atoms outside the condensate. This is an excellent 
approximation for a dilute Bose gas at low temperature 
where the depletion of the condensate is negligible. 

An important question to discuss concerning the 
ground state of a trapped Bose gas is the role of the 
interatomic potential. At first sight one would in fact 
expect that the role of interactions be negligible in a 
dilute system, where the usual expansion parameter a ^p 
is extremely small. Actually it turns out that the interac- 
tion can have a deep influence on the solution of Eq. (1), 
where its effect turns out to be fixed by the adime n- 
sional parameter Nala^o, where a^o = \/^/{m(OHo) is 
the harmonic oscillator length. This parameter can be 
indeed rather large despite the smallness of a^p. The 
final result is that the system is still fully Bose con- 
densed, but the structure of its wave function can be 
strongly affected by the interatomic forces. 

For the above reasons it is useful to discuss the solu- 
tion of the Gross-Pitaevskii equation in two relevant 
limits: the noninteracting model and the strongly repul- 
sive limit, which corresponds to the Thomas-Fermi ap- 
proximation. We will make a comparison with the exact 
numerical solution of the Gross-Pitaevskii equation in 
order to point out the role of the interaction. We will 
devote special attention to the structure of the conden- 
sate wave function near the boundary, close to the classi- 
cal turning point, and we will finally study the case of 
quantized vortices. 



2. The Noninteracting Model 

When the scattering length a vanishes, the problem 
reduces to the solution of a one-body Schrodinger equa- 



tion. Let us put a = in Eq. (1) and take the external 
potential as an anisotropic harmonic oscillator: 



KxtCr) = ^ ((ol x^ -\-a)ly^-\- (d\ z^). 



The ground state wave function becomes 



(2) 






^y.y^.xey 



,(3) 



where A = (ojco^. The Gaussian has different transverse 
and vertical widths. In particular one has 
(x") = (y') = {\l2)al and <z') = {\l2)k-'al. The chemi- 
cal potential is (1 -i- Xl7)fi(x>^ and coincides with the 
energy per particle, while the kinetic energy per particle 
has the simple form 



f = i<PV^(l.A/2)^«.. 



1 

2' 



(4) 



An interesting quantity to discuss is the ratio V(/?f )/(/?!) 
which provides a measure of the aspect ratio , character- 
izing the anisotropy of the velocity distribution. Using 
the wave function Eq. (3) one finds 



V(p,^)/(p.^) = V0OA?) = Va. 



(5) 



Values of the aspect ratio different from 1 reflect a 
peculiar and unique feature of Bose-Einstein condensa- 
tion. 



3. The Strongly Repulsive Limit (Thomas- 
Fermi Approximation) 

The opposite limit is obtained when the interaction is 
so strong, or the number of particles so large, that the 
kinetic energy term V^t/^ can be neglected in the Gross- 
Pitaevskii equation, Eq. (1). It corresponds to very large 
values of the dimensionless parameter 



Zi^aN 



(6) 



Also in this case the solution of the Gross-Pitaevskii 
equation is trivial and the wave function has the form: 

p(r) = lijjirr = 4^ [/^ - VUr)] (7) 

if the right hand side is positive, and p = elsewhere. 
The chemical potential is easily calculated by imposing 
the normalization condition J p 6r = N. One finds 
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IX- 



1 /15 \^^ 



(8) 



Due to the different scaling properties of the wave 
function with respect to th e variable z (compare Eqs. (3) 
and (7)), the aspect ratio ^ipDlipl) in this case is equal 
to A differently from the noninteracting case Eq. (5). 

The wave function Eq. (7) is expected to approximate 
well the exact solution of the Gross-Pitaevskii equation 
Eq. (1) for large A^, apart from the structure of the 
surface region where the exact wave function has to 
vanish smoothly. Some relevant observables, as the ki- 
netic energy, can be significantly affected by this sur- 
face structure, as we will see in Sec. 5. 



4. Solution of the Gross-Pitaevskii 
Equation 

The Gross-Pitaevskii equation, Eq. (1) can be solved 
numerically [5,6]. We transform the differential equa- 
tion in a functional minimization and use a steepest 
descent method to solve the minimization problem on a 
grid of points. As an example of atoms with repulsive 
interaction we choose ^^Rb, as in the experiment of Ref . 
[1]. For the 5 -wave triplet-spin scattering length we use 
a = 100^0, where ao is the Bohr radius. The asymmetry 
parameter is taken A = cojco^ = Vs and the axial fre- 
quency (oJ2tt = 220 Hz. The corresponding character- 
istic length is a^ = 1.222 X 10""^ cm. 

Results for the chemical potential and the energy per 
particle are shown in Table 1. Both quantities are ex- 
pressed in units of ^(o±, or of the equivalent temperature 
^(oJkB = 3.73 nK. The partial contributions to the en- 
ergy per particle coming from the kinetic energy (kin), 
the harmonic oscillator potential (HO) and the internal 
potential energy (pot) are also given. The A^ = 1 case 
coincides with the noninteracting anisotropic harmonic 

Table 1. Results for the ground state of ^^Rb atoms in a trap with 
(joJIti = 220 Hz and A = a)Ja)±_ = V8. Chemical potential and energy 
in units ^a)±_ 



N 


f^ 


(E/N) 


{E/N)un 


(E/N)uo 


(£W)pot 


1 


2.414 


2.414 


1.207 


1.207 


0.000 


100 


2.88 


2.66 


1.06 


1.39 


0.21 


200 


3.21 


2.86 


0.98 


1.52 


0.36 


500 


3.94 


3.30 


0.86 


1.81 


0.63 


1000 


4.77 


3.84 


0.76 


2.15 


0.93 


2000 


5.93 


4.61 


0.66 


2.64 


1.32 


5000 


8.14 


6.12 


0.54 


3.57 


2.02 


10000 


10.5 


7.76 


0.45 


4.57 


2.74 


15000 


12.2 


8.98 


0.41 


5.31 


3.26 


20000 


13.7 


9.98 


0.38 


5.91 


3.68 



oscillator: in this case the total energy per particle coin- 
cides with the analytic value (1 -H A/2) = 2.414. When A^ 
increases the repulsion among atoms tends to lower the 
central density, expanding the cloud of atoms towards 
regions where the trapping potential is higher. A typical 
profile of the condensate wave function ijj is plotted 
along the x-axis for A^ = 5000 in Fig. 1. The exact min- 
imization of the Gross-Pitaevskii functional (solid line) 
is compared with the noninteracting case (dashed line) 
and the Thomas-Fermi limit (dot-dashed). 
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Fig. 1. Ground state wave function (in arbitrary units) along x (in units 
aA_) for 5000 atoms of ^^Rb. Dashed line: noninteracting case. Dot- 
dashed line: Thomas-Fermi approximation. Solid line: numerical so- 
lution of the Gross-Pitaevskii equation. 

When A^ is large we observe an increase of both 
interaction and harmonic oscillator potential energy per 
particle (the latter effect follows from the expansion of 
the cloud). Conversely, the kinetic energy per particle 
decreases because the density distribution is flattened. 
In the strongly repulsive limit, N ^ ^, one should find 
that the internal potential energy is much greater than 
the kinetic energy. Indeed the convergence towards this 
limit turns out to be rather slow as we will show in the 
next section. 

Another interesting quantity which can be easily cal- 
culated from the ground state wave function is the aspect 
ratio of the velocity distribution, that is the ratio ^{pl)l 
ipl). This quantity is equal to Va in the noninteracting 
case and should approach A in the strongly repulsive 
limit. The numerical results, as a function of A^, are 
shown in Fig. 2. The two limiting cases are shown as 
dashed lines. One clearly sees that the convergence to 
the value 2.828 = A is very slow; the aspect ratio remains 
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Fig. 2. Ratio of the axial to transverse average velocity as a function 
of A/^in ^^Rb. The lower and upper dashed lines corresponds to V A and 
A, respectively 



■^ 0.3 - 




well below the asymptotic value even for A^ = 20000. 
The aspect ratio measured in Ref. [1] is estimated to be 
about 50 % larger than the noninteracting value, while 
the number of particles is of the order of 5000. The 
agreement with our results is good, even if one has to 
consider that the experimental estimate implicitly as- 
sumes a ballistic expansion of the atoms after switching 
off the external trap. The effects of the interaction on the 
expansion of the gas should be explicitly taken into 
account in order to draw more definitive conclusions. 

As an example of atoms with attractive interaction we 
choose ^Li, as in the experiment of Ref. [2]. For the 
5 -wave triplet-spin scattering length we use a= — 27^0. 
The axial frequency reported in Ref. [2] is cjjln =117 
Hz and the corresponding characteristic length is 
ax_ = 2.972 X 10""^ cm. The transverse frequency is oj^/ 
217=163 Hz, so that the asymmetry parameter is 
A = (jjjoix^ = 0.72. 

The first important point to stress is that Gross- 
Pitaevskii functional has no global minimum for nega- 
tive scattering length. This reflects the tendency of the 
system to collapse. For spatially inhomogeneous sys- 
tems, however, the zero-point energy can exceed the 
attractive potential, producing local minima of the func- 
tional when the density of atoms is not too high. 

The most striking difference with respect to the re- 
pulsive case is that here the central density of the cloud 
increases rapidly with A^, as shown in Fig. 3. This is the 
effect of adding more and more attractive potential en- 
ergy. When the central density reaches a certain critical 
limit the system collapses and the solution of the Gross- 
Pitaevskii equation does not converge anymore. In ^Li, 
with the input parameters given above, the critical num- 
ber A^ turns out to be about 1400. In [6] we have found 



Fig. 3. Ground state wave function (in arbitrary units) along x (in units 
a^} for 1000 atoms of ^Li. Dashed line: noninteracting case. Solid line: 
numerical solution of the Gross-Pitaevskii equation. 

that a stationary solution of the Gross-Pitaevskii equa- 
tion with larger values of A^ can be obtained if a vortex 
line is present in the system. The possible occurrence of 
vortices in these trapped Bose gases will be discussed in 
the Sec. 6. 



5. Wave Function at the Boundary 

As one clearly sees in Fig. 1 , the Thomas-Fermi ap- 
proximation fails to reproduce the structure of the order 
parameter at the surface of the atomic cloud in the case 
of positive scattering length (repulsive interaction). Sev- 
eral measurable quantities can be significantly affected 
by the behavior of the wave function in this region. In 
order to provide a good model for these quantities one 
has to go beyond the Thomas-Fermi approximation. One 
of these relevant observables is the kinetic energy 



'"2mJ 



drIVi/rll 



(9) 



In fact the Thomas-Fermi approximation (7) for the 
wave function is not appropriate to evaluate £kin; it pro- 
duces a logarithmic divergence in the integrand of Eq. 
(9), occurring at the classical turning point, where 
Ve^^ = fx . This reveals that the evaluation of £kin requires 
higher accuracy in the description of the boundary re- 
gion. In order to provide the proper description of 
the condensate wave function near the boundary we 
have recently proposed [7] a suitable expansion of the 
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Gross-Pitaevskii equation near the classical turning 
point. The resulting analysis allows one to obtain the 
proper expansion for the kinetic energy. We briefly 
sketch here the case of isotropic traps (a^ = «z = «ho)- 
Let R be the boundary of system spherical system, 
determined by the equation fx = V^^tiR)- Near this point 
one can carry out the expansion 



VUr) - fj.^(r-R)F 



(10) 



where F = moj^oR is the modulus of the attractive exter- 
nal force evaluated at r = R. Close to the boundary, 
where \r — R\ <kR, the Gross-Pitaevskii equation takes 
the form 






R)Fifj + 



2m 



il/^ = 0. (11) 



solution of the Gross-Pitaevskii equation in the external 
surface profile. An example is given in Fig. 4 for 10^ 
atoms of ^^Rb in a spherical trap. 



0.06 



0.04 



-^ 



0.02 - 



Let us now introduce the dimensionless variable 




(r-R) 



(12) 



where 



H> 



(13) 



is a typical thickness of the boundary giving, as we will 
see later, the distance from the classical radius R where 
the Thomas-Fermi approximation starts failing. Then 
we introduce the dimensionless function cf) defined by 



^(r) = 



1 



diSiraY' 



^(a 



(14) 



in terms of which the Gross-Pitaevskii equation, Eq. 
(11) takes the universal form 



r -{^^ <!>')<!> = 0. 



(15) 



Its solution provides, via Eqs. (12-14), the proper struc- 
ture of the condensate wave function near the classical 
turning point 7? . It is worth noting that Eq. ( 1 5) does not 
depend on the form of the external potential nor on the 
size of the interatomic force. These physical parameters 
enter the transformations Eqs. (12) and (14) which fix, 
together with the solution of Eq. (15), the actual behav- 
ior of the wave function ijj. Equation (15) can be solved 
numerically. The function cf) behaves like v— ^ for 
^ -> - cx) and like T''^ exp[- (2/3)^'''] in the opposite 
limit ^ ^ ^. The corresponding condensate wave func- 
tion ijj matches the Thomas-Fermi wave function at the 
left of the classical turning point and follows closely the 



Fig. 4. Condensate wave function (arbitrary units) for 10^ atoms of 
^^Rb in a spherical harmonic trap of length ^ho- Solid line; numerical 
solution of the Gross-Pitaevskii equation, Eq. (1). Dot-dashed line: 
Thomas-Fermi approximation, Eq. (7) (indistinguishable from the 
solid line in the inner part). Dashed line: surface profile obtained from 
the universal equation, Eq. (15). 

The kinetic energy can be calculated by matching 
properly the Thomas-Fermi approximation and the solu- 
tion of the universal equation, Eq. (15). This yields the 
result 

where the radius R is related to N by the equation 
R' 



N = 



15 a auo ' 



(17) 



Equation ( 1 6) provides the proper behavior of the ki- 
netic energy in the large N limit where R » Uho- In Fig. 
5 we compare the results obtained from Eq. (16) and 
from the Gross-Pitaevskii equation. One can see that the 
convergence is reached for relatively large values of N. 



6. Vortices 

The structure of vortices in a trapped Bose gas can be 
naturally investigated in the present formalism. Let us 
consider states having a vortex line along the z -axis and 
all the atoms flowing around it with quantized circula- 
tion. One can write the axially symmetric condensate 
wave function in the form 
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0.3 




Fig. 5. Kinetic energy per particle, in units ^(Who, for ^^Rb in a 
spherical harmonic trap as a function of the number of condensed 
atoms. Solid line: from the solution of the Gross-Pitaevskii equation 
(1). Dashed line: approximation Eq. (16). 



^(r) = ij/{r) cxp[iS{r)] 



(18) 



where il/(r) = vp(r) is the modulus, while the phase S 
acts as a velocity potential: v = (^/m)VS. By choosing 
S = Kcj), where (p is the angle around the z-axis and k is 
an integer, one has vortex states with tangential velocity 



h 

V= K, 

mr± 



(19) 



with rl=x^ -\- y^. The number k is the quantum of circu- 
lation, and the angular momentum along z is Nk^ . 

If the complex wave function Eq. (18) is used in the 
derivation of the Gross-Pitaevkii equation, one gets 






477^ a . , . . ,2 

+ \ih(r)r 

m 



il/(r) = fxil/(rX 



(20) 



which differs from Eq. (1) only for the addition of a 
centrifugal potential. This new term forces the solution 
ijj to vanish on the z-axis for k i^ 0. 

For noninteracting particles one falls again in the case 
of the stationary Schrodinger equation for the an- 
isotropic harmonic potential. For instance the /<: = 1 so- 
lution has the form 



il/(r) cc n exp 



il^^'-'^'^ 



(21) 



To get the energy per particle for the k i^ states one 



has simply to sum k^o)± to the energy per particle of the 
ground state without vortices. 

In the interacting case the kinetic energy can not be 
neglected even for large N, since it determines the struc- 
ture of the vortex core. In particular, the balance be- 
tween the kinetic energy and the interaction energy fixes 
a typical distance over which the condensate wave func- 
tion can heal. For a dilute Bose gas the healing length is 
given by 



^=(87rpa)- 



(22) 



where p is the density of the system. In the case of a 
vortex it corresponds to the distance over which the wave 
function increases from zero, on the vortex axis, to the 
bulk density. For the trapped atoms in the A/^ — > co limit 
one finds 



R \r) 



(23) 



Thus the healing length is small compared with the size 
of the cloud if R is much bigger than a^. 

The critical angular velocity required to produce vor- 
tex states is easily calculated once the energies of the 
states with and without vortices is known. One has to 
compare the energy of a vortex state in frame rotating 
with angular frequency /2, that is (E — flL^), with the 
energy of the ground state with no vortices. Since the 
angular momentum per particle is k/i, the critical angu- 
lar velocity is 



n,= {M)-\{EIN\- {EIN\l 



(24) 



In the noninteracting case the difference of energy per 
particle is simply Kho}^, so that /2c = w^. 

We have solved numerically the Gross-Pitaevskii 
equation (20) both for rubidium and lithium. In Fig. 6 
we show the wave function of a cloud of 5000 ^^Rb 
atoms; the /<: = 1 wave function (Fig. 6b), which corre- 
sponds to atoms flowing around the z-axis with angular 
momentum Nh , is compared with the /<: = ground state 
(Fig. 6a). The atoms are pushed away from the axis 
forming a toroidal cloud. From the energy of the vortex 
states we calculate the critical angular velocity, through 
Eq. (24). The results for /<: = 1 are shown in Fig. 7. The 
critical angular velocity decreases rapidly with N. For 
N > 5000 it is less than 40 % of the noninteracting value, 
given by the transverse angular frequency Wi of the 
trap. The healing length is the distance over which the 
wave function grows from zero to the bulk value. In the 
limit of large systems it can be approximated by Eq. (22) 
with p equal to the density in the central part of the 
toroidal distribution. Both the estimate of ^ and fl^ 
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Fig. 6. Wave function, in arbitrary units, of 5000 Rb atoms. Spatial 
coordinates in units a^. a) Ground state, b) Vortex state with k=\. 

obtained in this way are in qualitative agreement with 
the behavior of the numerical solutions. 

Coming back to the question of the stability for nega- 
tive scattering length, we notice that, when the local 
minimum associated with wave functions of the form 
shown in Fig. 3 disappears, nothing prevents a /7n6>n the 
existence of other local minima associated with different 
configurations. Such configurations should have local 
density lower than the critical one. A natural way to 
obtain a favourable situation is to move the atoms away 
from the z -axis, conserving the total number of particles. 
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Fig. 7. Critical angular velocity, in units 6>i, for the formation of 
K = 1 vortices in ^^Rb vapor as a function of A^. 

This happens in the presence of a vortex. In Fig. 8 we 
show the wave function for 1000 ^Li atoms with no 
vortices (Fig. 8a) and with an axial vortex of unit circu- 
lation (Fig. 8b). We use the same units in both cases, so 
one can see that the maximum value of the wave func- 
tion inside the toroidal distribution of the vortex is ap- 
proximately a factor two lower than the central value in 
the state with no vorticity (the density is four times 
smaller). The critical angular frequency for the forma- 
tion of the vortex state in Fig. 8 is 1.12 times the trans- 
verse angular frequency of the trap. In systems with 
attractive interaction the critical angular velocity is 
larger than for noninteracting particles, while the oppo- 
site is true for repulsive interaction. This is because it 
costs internal potential energy to lower the average den- 
sity, as the vortex does, for attractive interactions. How- 
ever, once a vortex is created, the corresponding state is 
more stable than in the absence of vorticity: one can put 
more atoms inside the rotating cloud before reaching the 
critical density for the final collapse. Indeed we find 
local minima of the Gross-Pitaevskii functional for A^ 
much larger than 1400 if /c> 0. For /c = 1 we find a 
critical value of N — 4000; for /c = 2 and 3 we find 
critical values of 6500 and 8300, respectively. It is worth 
mentioning that the number of particles in the conden- 
sate reported in the experimental work of Ref . [2] is an 
order of magnitude higher than the critical value for the 
stability of the Gross-Pitaevskii solution without vortic- 
ity (A^ — 1400). The presence of vortices might explain 
the large size of the observed Bose-condensed gas. Fur- 
ther experimental data are needed to draw more defini- 
tive conclusions. 



543 



Volume 101, Number 4, July-August 1996 

Journal of Research of the National Institute of Standards and Technology 



About the authors: Franco Dalfovo is a member of the 
research staff at the Department of Physics, University 
of Trento, Italy. Lev Pitaevskii is a Visiting Professor at 
the Department of Physics of the Technion, Haifa, Israel, 
and a Main Researcher at the Kapitza Institute for Phys- 
ical Problems, Moscow, Russia. Sandro Stringari is a 
professor of theoretical physics at the Department of 
Physics, University of Trento, Italy. 




2 

X 



Fig. 8. Wave function, in arbitrary units, of 1000 ^Li atoms. Spatial 
coordinates in units a±. a) Ground state, b) Vortex state with k= I. 
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